#!/usr/bin/perl
# 2012/11/16 Checking fastq quality (modified from priyankara version) Takes long time.

use strict;
use warnings;
use IO::File;
use Getopt::Std;
use Carp qw(croak);
$| = 0;

if(@ARGV != 4 or $ARGV[0] eq '--help'){
	die "Usage : $0 -i [fastq file] -o [output pdf file]\n";
}

#---------------------------------------
# parameter
#---------------------------------------
my %opt;
getopts("i:o:", \%opt);
my $FILE_fastq = $opt{i};
my $FILE_pdf = $opt{o};

my $FILE_tmp = 'tmp_' . $FILE_fastq;

#---------------------------------------
# read quality values
#---------------------------------------
my %data;
{
	my $fh_in = IO::File->new($FILE_fastq) or die "cannot open $FILE_fastq: $!";

	my $line = 0;
	my $quality;
	my $sequenceNum = 0;
	while($_ = $fh_in->getline()){
		s/\r?\n//;
		$line++;
		if($line % 4 == 1){
			$sequenceNum++;
			if($sequenceNum % 5000000 == 0){
				&Progress("$sequenceNum finished...");
			}
		}elsif($line % 4 == 2){
		}elsif($line % 4 == 0){
			my @characters = split //;
			for(my $i = 0; $i < @characters; $i++){
				my $Q = ord($characters[$i]) - 33;	# quality value
				$data{$i}{$Q}++;
			}
		}
	}
	$fh_in->close();
}


#---------------------------------------
# output tmp file
#---------------------------------------
open OUT, ">$FILE_tmp" or die "cannot write $FILE_tmp: $!";
foreach my $k(sort {$a <=> $b} keys %data){
	my $flag = 0;
	print OUT "p$k <- c(";
	foreach my $Q(keys %{$data{$k}}){
		if($flag == 0){
			if($data{$k}{$Q} == 1){
				print OUT "$Q";
			}else{
				printf OUT "rep(%d,%d)", $Q, $data{$k}{$Q};
			}
			$flag = 1;
		}else{
			if($data{$k}{$Q} == 1){
				print OUT ", $Q";
			}else{
				printf OUT ", rep(%d,%d)", $Q, $data{$k}{$Q};
			}
		}
	}
	print OUT ")\n";
}
print OUT "pdf(\"$FILE_pdf\")\n";
print OUT "boxplot(";
foreach my $k(sort {$a <=> $b} keys %data){
	print OUT "p$k, ";
}
print OUT<<"EOF";
xlab="read positions", ylab="Quality score", main="$FILE_fastq")
dev.off()
EOF

close OUT;

system("Rscript --vanilla --slave $FILE_tmp") == 0 or die "R failed: $!";
unlink($FILE_tmp);


sub Progress{
	my ($message) = @_;
	my $time = &time();
	warn "[$time] $message\n";
}


sub time{
	my ($sec, $min, $hour, $mday, $mon, $year) = localtime(time);
	my $time = sprintf ("%d/%0d %02d:%02d",$mon+1, $mday, $hour, $min);
	return $time;
}

